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Abstract. An adhesive unilateral contact of elastic bodies with a small viscosity in the 
linear Kelvin-Voigt rheology at small strains is scrutinized. The flow -rule for debonding the 
adhesive is considered rate -independent and unidirectional, and inertia is neglected. The 
asymptotics for the viscosity approaching zero towards purely elastic material involves a 
certain defect-like measure recording in some sense natural additional energy dissipated 
in the bulk due to ( vanishing) viscosity, which is demonstrated on particular 2-dimensional 
computational simulations based on a semi-implicit time discretisation and a spacial dis- 
cretisation implemented by boundary-element method. 



1. INTRODUCTION, QUASISTATIC DELAMINATION PROBLEM 

Quasistatic inelastic processes on surfaces of (or interfaces between) solid elastic bodies 
like fracture or delamination (or debonding) of adhesive contacts have received intensive 
engineering and mathematical scrutiny during past decades. Often, the time scale of such 
processes is much faster than the external loading time scale, and such processes are then 
modelled as rate independent, which may bring theoretical and computational advantages. 
Yet, the above mentioned inelastic phenomena typically lead to sudden jumps during evo- 
lution, which is related with the attribute of nonconvexity of the governing stored energy 
(cf. here the non-convex term J r |zKm m dS in ( l6dl i below), and then it is not entirely clear 
which concept of solutions suits well for the desired specific application. 

The "physically" safe way to coup with this problem is to reduce rate-independency 
on only such variables with respect to which the stored energy is convex, the resting ones 
being subjected to certain viscosity (or possibly also inertia). Here we neglect inertia from 
the beginning, which is addressed as a quasistatic problem; cf. ll29[ Sect. 5] for the dy- 
namical case. Moreover, such viscosities can be (and in most materials also are) very 
small, and in engineering literature are almost always completely neglected. However, al- 
though arbitrarily small, such viscosities are critically important to keep energetics valid. 
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It therefore makes a sense to investigate the asymptotics towards purely elastic materials 
when these viscosities vanish. In the limit, we thus get some solutions of the underlying 
rate-independent system which, however, might (and, in specific applications, intentionally 
should) be different from solutions arising when viscosity are directly zero and global- 
energy-minimization principle is in play, cf. also Remark|2]below. 

In this article, we will confine ourselves to visco-elastic bodies at small strains and 
we consider the viscosity in the Kelvin-Voigt rheology, which is the simplest rheology 
which makes the desired effect of natural prevention of the too-early delamination, cf. 
1 2911 . It should also be emphasized that our viscosity is in the bulk while the inelastic 
delamination itself is considered fully rate-independent, in contrast to a usual vanishing- 
viscosity approach as e.g. in [0, El , 16. 21, 311. A certain bulk viscosity but acting 
on displacement itself rather than on the strain was considered in J2]. 

For notational simplicity, we consider a single visco-elastic body occupying a bounded 
Lipschitz domain Q, c M. d and the adhesive contact on a part T c of the boundary d£l, so that 
we consider dQ. = F c U F D U T N U with disjoint relatively open T c , T D , and T N subsets of 
(9Q and with N having a zero (c/-l)-dimensional measure. All results are, however, valid 
equally for delamination on boundaries inside O, i.e. an adhesive contact between several 
visco-elastic bodies. For readers' convenience, let us summarize the notation used below: 



d dimension of the problem id = 2, 3), 
if displacement (defined on fl), 
z delamination parameter (defined on T c ), 
e(u) = 4(Vw) T + jVu small-strain tensor, 
C tensor of elastic moduli of the 4 th -order, 
X a "Kelvin-Voigt" relaxation time, 
xC viscous-moduli tensor, 
e = eixu + u) 

K the matrix of elastic moduli of the adhesive, 



E Young modulus, 

y Poisson ratio, 

a fracture toughness, 

driving energy for delamination, 

t traction stress vector (acting on 

r N u r c ), 

t„, t, normal or tangential compo- 
nent of t, 

/ bulk load (acting on CI), 

g surface load (acting on T N ), 

w D surface displacement loading 

(on r D ). 



Table 1 . Summary of the basic notation used thorough the paper. 



We consider the standard model of a unilateral frictionless Signorini contact. The qua- 
sistatic boundary- value problem for the displacement u on Q. and the so-called delamination 
parameter z on F c valued in [0, 1], representing Fremond's concept [1JJ| of delamination, 
considered in this paper is: 



with 



e = e(u, ii) — x e (u) + e(u) 



divCe + / = 
u - w D 

m = 8 

t t (e) + z(Ku-((Ku)-n)n) = 0, 

u-ft > 0, t n (e)+z(Ku)-n > 0, (t n (e)+z(Ku)-n)(u-n) = 0, 
z < 0, d < a, z(d-a) = 0, 

3 e \Ku-u + AfoiKz) 



on il, 
on F D , 
on F N , 



on r r 



(la) 
(lb) 
(lc) 

(Id) 



where we use the usual "dot-notation" for the time derivative, i.e. (■)'=§}, where the set- 
valued mapping A^[o,ij : M =4 K assigns z e R the normal cone A^[o,i](z) to the convex set 
[0, 1] c R at z e M, and where the traction stress and its normal and tangential components 
are defined on F c U F N respectively by the formulas 



Ke) = (Ce)\ r H, t n (e) = ti(Ce)\ r ft, t t (e) = (Ce)| r # - t n (e), (le) 



where it = n(x) is the unit outward normal to F := c?Q. We further consider the initial-value 
problem for (QJ-e) by prescribing the initial conditions 



w(0) = m and z(0) = z Q . (If) 

Of course, the loading /, g, and w D in (TjJ-c) depend on time t. The parameter a > in 
( TTdl i is a given phenomenological number quantity (possibly as a function of x e T c ) with a 
physical dimension J/m 2 with the meaning of a specific energy needed (and thus deposited 
in the newly created surface) to delaminate lm 2 of the surface under adhesion or, equally, 
the energy dissipated by this delamination process; in fact, ([8]) below reflects the latter 
interpretation. In engineering, a is also called fracture toughness (or fracture energy). 

As already mentioned, in conventional materials, the viscosity and thus relaxation time 
X > is mostly very small in comparison with external force loading time-scale, and it is 
worth studying the asymptotics forx — > 0. Formally, the inviscid limit problem arising for 
X — » is a quasistatic problem for purely elastic material which consists in replacing ([Tal l 
by 

div Ce(u) + f = on Q, (2) 

and in replacing t(e) by t(e(u)) in dTcb and similarly t n (e) and t t (e) by t n (e(u)) and t t {e(u)) 
in dTdb with t(-), t n (-), and t t (-) again from ([Tel l. This limit rate- independent problem itself, 
however, does not record any trace of energy dissipated by viscosity in the bulk during 
rupture of the delaminating surface, but there are explicit examples, cf. l29ll . showing that 
this energy is not negligible no matter how the viscosity coefficient ;f > is small, which 
leads to a notion of Kelvin-Voigt approximable solution to this limit rate-independent prob- 
lem involving a certain, so-called defect measure recording the "memory" of this dissipated 
energy which somehow remains even if viscosity coefficient ;f vanishes (i.e. is passed to 0). 

The plan of the paper is as follows: First, in Section [2] we formulate the above initial- 
boundary-value problem ([T} weakly and briefly present the main results about a-priori es- 
timates and convergence for x — > to the inviscid quasistatic rate-independent problem, 
leading to the above mentioned approximable solutions and defect measures, mainly taken 
from I29I1 . In Sect. [3] we perform time discretisation by a semi-implicit scheme and present 



some convergence results again from H29H . and prove that the residuum in the discrete en- 
ergy balance converges to zero if the time step goes to 0. Merging Sections [2] and [3] 
this energy-residuum convergence serves as an important ingredient for controlling conver- 
gence of the discretisation with dependence on convergence of viscosity. This is eventually 
used in Section|4]where, making still a spacial discretisation by boundary-element method 
(BEM), we perform computational experiments both with a one -dimensional example from 
02911 with a known solution to tune parameters of the algorithm and eventually with a non- 
trivial two-dimensional example. In this last example, we (to our best knowledge histori- 
cally for the first time) present numerical study of a nontrivial, spatially non-homogeneous 
defect measure. 

2. INVISCID PROBLEM AS A VANISHING-VISCOSITY LIMIT 

The weak formulation of the initial-boundary value problem ([TJ is a bit delicate due to 
the doubly-nonlinear structure of the flow rule for z on F c without any compactness (i.e. 
without any gradient theory for z) and with both involved nonlinearities unbounded due to 
the constraints z < and z > (while the third constraint z < 1 is essentially redundant if 
the initial condition satisfies it). This would make serious difficulties in proving the exis- 
tence of conventional weak solutions. Benefiting from rate-independency of the evolution 
rule for z, we can cast a suitable definition by combining the conventional weak solution 



concept for u and the so-called energetic-solution concept II 18112211 of z as in 12711 . 

Considering a fixed time horizon T > 0, we use the shorthand notation / = (0, T), 
I = [0,71, Q = Ix£l, Q = 7xQ with £1 the closure of £1, 2 D = 7xF D , andS c =/xf c . We 



will assume, without substantial restriction of generality of geometry of the problem, that 

dist(F c , r D ) > 0, meas^i(r D ) > 0, meas,/_i(r c ) > 0. (3) 

We first make a transformation of the problem to get time constant Dirichlet condition. 
To this goal, we first consider a suitable prolongation m d of w D defined on Q, i.e. m d |e d = vv D . 
Then we shift u to u + u D , and rewrite ([TJ for such a shifted u. Thanks to the first condition 
in (01, we can assume that m d |e c = so that ( [Tdb remains unchanged under this shift. The 
equations (QJ-c) transform in such a way that the original loading /, g, and w D as well as 
the original initial data uq are respectively modified as follows: 

/ replaced by / + divCe D with e D = e(xu D +u D ), (4a) 

g replaced by g + (Ce D )|r N «, (4b) 

w D replaced by 0, (4c) 

uq replaced by uq - u D (Q). (4d) 

We will use the standard notation W l ' p (Sl) for the Sobolev space of functions having 
the gradient in L p (Q.;R d ). If valued in W with n > 2, we will write W l - P (Q:,W), and 
furthermore we use the shorthand notation H l (£l; W 1 ) = W l 2 (Cl; R"). We also use the no- 
tation of " ■ " and " : " for a scalar product of vectors and 2nd-order tensors, respectively. 
Later, Meas((5) = C(Q)* will denote the space of measures on the compact set Q. For a 
Banach space X, L P (I;X) will denote the Bochner space of X-valued Bochner measurable 
functions u : I — > X with its norm ||m(-)II in L P (I), here || ■ || stands for the norm in X. 
Further, BV(I; X) will denote the space of mappings u : I — » X with a bounded variations, 
i.e. supQ^^,^ <t l<t <t YIi=\ \\ u (U)—u(ti-\)\\ < oo where the supremum is taken over all fi- 
nite partitions of the interval [0, T]. Also, we will use H X {I;X) for the Sobolev space of 
X-valued functions with distributional derivatives in L 2 (I; X). To accommodate the trans- 
formation (0]i into the weak formulation, we introduce the functional f(f) € H l (Q;M. d )* 
by 

<f(/),v):= f f{t)-v-Ce(xu D {t)+u {t)):e{v)+ f g(t)-vdS. (5) 
Jn Jr N 

Definition 1. The couple (u x ,z x ) with u x e H l (I;H l {Q.;M. d )) and z x 6 BV(I;L\T C )) n 
L°°(S C ) is called an energetic solution to the initial-boundary-value problem (Q~|i under the 
transformation (0} if 

(i) the momentum equilibrium (together with Signoring boundary conditions) in the weak 
form r r-T 

I Ce(xu x +u x ):e(v-u x ) dxdt + I z x Ku x -(v-u x ) dSdt > I (j(t),v-u x ) dt (6a) 
Jq J^c Jo 

with f defined in © holds for any v e H l (I; H l (Q\ R d )) with v| Sc -« > and v|s D = 0, 

(ii) the so-called semi-stability of the delamination holds for any t e [0, T]: 

Ku x (t, x)-u x (t, x) < 2a(x) or z x (t,x) - for a.a. x e F c , (6b) 

(iii) and the energy equality 

&(t,u x (t),z x (t))+ f ^ x^ e (u x ):e(u x ) dxdt + f a(z -z x (t))dS = £{0,uq,zo) + f (f,u x )dt 
Jo Jn Jr c Jo 

(6c) 

holds for any t 6 [0, T] with f defined again by (0, and with 

f -Ce(u):e(u)dx - <f(f), «> + f -zKu-udS if u-n > 0, < z < 1 on F c , 
Jn2 Jr c 2 

and if u = on r D , 
+oo else, 

(6d) 



£{t, u,z) : 



(iv) the initial conditions (fTft understood transformed as d4cfl i hold. 

This definition is indeed well selective in the sense that any smooth energetic solution 
solves also (Q]i in the classical sense. Due to d6cT i, ${t,u x {t),z x {t)) < oo so that it holds 
u x \z. c -n > 0, u x \j, D = 0, and < z x < 1 if the initial conditions satisfies these constraints 
so that <?(0, U(),zo) < °°. Note also that ([T} has an abstract structure of the initial-value 
problem for the triply nonlinear system of two evolution inclusions: 

\8% x \u + d u S{t, u,z)3 0, m(0) = M , (7a) 
d- z M x {'z) + d-S(t, u, z) 3 , z(0) = zo, (7b) 

with [■]' denoting the (partial) Gateaux differentials and d denoting the partial subdif- 
ferentials in the sense of convex analysis, with S from (16dt . and with the ^--dependent 
(pseudo)potential of dissipative forces 3t x given by 



Au,z) 



I — Ce(u):e(u) dx + j a\z\dS if z < a.e. on T c , 
Ja 2 Jr c (8) 

+oo otherwise. 



Also note that d6bb is equivalent to the integrated form of the abstract semistability S (t, u x (t), z x (t)) < 
<o(t,u x (t),z) + &Q(z-z x (t)) holding for any z > 0, where we wrote briefly Mo(u,z) = 
& x (0, z) =: &o{z). This means here: 

V2eL°°(r c ), 0<2<z^(f): { {z x (t)-z){Ku x (t)-u x (t) -2a)AS <Q. (9) 

Jr c 

We will generally assume the following data qualification: 

C > (= positive definiteness), (10a) 

f e W 1A (r,H\Cl;R d )*), u Q e H\Cl;R d ), z e L°°(r c ), (10b) 

uo\r c -ft > 0, < zo ^ 1 a.e. on F c , and (10c) 

K.Uo(x)-uq(x) < 2a or zo(x) = for a.a. x € F c . (10d) 

Note that the qualification of f in dlObl ) represents, in fact, assumptions on /, g, and w D 
in the original problem ([TJ, and that (TTOb.d) means semi-stability of the initial condition 
(uq,zq). Under ([Tol l, existence of the solutions due to Definition Q] ca n, in fact, be proved 
by limiting the discrete solutions ( [Pil l, cf. Lemma[TJand details in |27, 29ll . 



Proposition 1 (Vanishing viscosity limit, 12911 ). Let (0 and ( ITOb hold, and let^ > 0. Then: 
(i) Any solution (u x ,z x ) according to Definition[T]satisfies the a-priori estimates: 



\"ALH,;HHnm)- C/ ^> (lla) 

l"^IL'»(/;//i(n;K'')) - C ' ^ llb ^ 

<C (11c) 



lr^llL»(Sc) n BV(/;L'(rc)) 

with C independent of^-. 

(ii) There are u e L°°(7; H \Q,; R d )), z e BV(I; L l (T c )), and /x eMeas(0, and a subsequence 
such that, forx ^ 0, 

u x (t)^>u(t) in// 1 (Q;R"') fora.a. te[Q,T], (12a) 

z x (t)^z(t) inL°°(r c ) for all te[0,T], (12b) 

xCe(u x ):e(u x ) — » /i inMeas(0. (12c) 

(iii) Any triple (u,z,fi) obtained by this way fulfills, for a.a. t 6 [0, T], the momentum 
equilibrium in the weak form, i.e. 

I Ce(u(t)):e(v-u(t))dx+ \ z(f)Ku(f)-(y-u(f)) dS > (f(t), v-«(/)> (13a) 
Jn Jr c 



for all v 6 H l (Q.; M. d ) with v\r c -ft > , furthermore the semi-stability 



Ku(t, x)-u(t, x) < 2a(x) or z(t, x)-0 fora.a. xeY c (13b) 
and eventually the energy equality 

<g(t,u{t),z{t)) + fa(z -z(t))dS + f [ji(dxdt) = S'(0,uo,zo)+ f <f,«)df. (13c) 
Jr c Jo Jn Jo 

The above assertion suggests the following: 

Definition 2. A triple (u,z,fi) with u e L°°(I;H l (£l;M. d )), z e BV(I\L l (T c )), and p e 
Meas(g), yu > 0, is called a Kelvin- Voigt-approximable solution to the quasistatic rate- 
independent delamination problem (Q~|i with x — transformed by © if (fT3l holds for 
a. a. t e /, and z(0) = Zo> an d if (u,z,fi) is attainable by a sequence of viscous solutions 
{(u x ,z x )} x >o in the sense (Tl2l i. 

The measure e Meas(0, invented in l29ll . occurring in Proposition [TJ represents a 
certain additional energy distributed over Q specified rather implicitly by (I12ct but anyhow 
with a certain physical justification. Similar concept has been invented in various other 
problems in continuum mechanics (particularly of fluids) under the name of defect mea- 
sures to reflect a possible additional energy dissipation of solutions lacking regularity and 
exhibiting various concentration effects in contrast to regular weak solutions where the de- 
fect measure vanishes, cf. era. Here, p. reflects the possible additional dissipated 
energy of Kelvin- Voigt-approximable solutions comparing to the so-called energetic solu- 
tions, cf. also Remark |2]below. A nonvanishing ji is vitally important and rather desirable 
in the context of fracture mechanics in contrast to the mentioned fluid-mechanical appli- 
cations where the phenomenon of nonvanishing ju is related "only" to a possible lack of 
regularity of weak solutions and is not entirely clear whether it has some physical justifica- 
tion and supported experimental evidence. 



3. TIME DISCRETISATION, CONVERGENCE 

Some solutions to the initial -boundary value problem (Q3 in accord to Definition Q] can 
be obtained rather constructively by a semi-implicit time discretisation. To facilitate the 
a-priori estimates, we again consider the transformation (HJi. Using an equidistant partition 
of the time interval [0, T] with a time step t > such that T/t € N, we consider: 



,k „k-l 



div Cef + f k T = with 4 = X e( r ~ ^ ) + e{u\) 

Llj — 

i(4) = 8r 

l t (e k T ) + z k - l (Ku k T -((K U k T )-fI) f T) = 0, 

u k T -n > 0, i n (e k T )+z k T - 1 (Ku k yri > 0, (t n (4)+4- 1 (K^)-n)(«r-n) = 0, 



z k < z k - x 



d k € 1 



t)* < a, 

k\ 



(4 - 4~ l m - a) = 0, 



2 K« + N m] (4) 



on £1, 
on T D , 

on r N , 



on r r 



(14a) 

(14b) 
(14c) 

(14d) 



with t(-), tnOX an d t t (-) from ([Tel l and with /* = f(kr) and g k = g(kr) with / and g from 
©, and proceeding recursively for k = 1, T/t with starting for k = 1 from 



U T = UQ 



and 



z?=zo. 



(15) 



The adjective "semi-implicit" is related with usage of z k ~ l in the first complementarity 
problem in ( I14db . instead of z k which would lead to a fully implicit formula. Such usage 
of z k ~ x leads to the decoupling of the problem: first we can solve (Tl4k-c) with the first 



complementarity problem in (I14db for u k and only after the rest of (I14dt for z k ; in fact, 
this can be understood as a fractional-step method, cf. also lEil Remark 8.25]. In addition, 
we can employ the variational structure of both decoupled problems. We thus obtain two 
convex minimization problems: first, we are to solve 

minimize g(kx, u, z k ~ l ) + tM x [ U ~ Ut ,o) 
subject to u e H l (Q.; R d ), u\r D = 0, u\ Tc -n > 

and, denoting its unique solution by u k , then we solve 

minimize S{kx, u k T , z) + M x (0, z-z 1 ^ 1 ) 
subject to zeL°°(r c ), < z < 4 _1 

with the stored energy § and the dissipation (pseudo)potential M x defined here by 

£{t,u,z)= f ^Ce(u):e(u)dx + f \z&u-udS - (f(t),u), (17a) 
Jn 2 Jr c 2 

M x {u,'z)= \ ?rCe(u):e(u)dx- \ a'zdS. (17b) 
Jn 2 Jr c 

Note that the constraints u\r c -n > 0, < z < 1, and i < 0, originally contained in $ and 
ffl x in (l6db and ||S), are now included in (TToT l so that we can equivalently use the smooth 
functionals S(t, •, •) and 0l x in ( fTTI i. Also note that M x (u, •) is degree-1 homogeneous so 
that the factor r does not show up in the functional in ( I16bb . in contrast to the degree-2 
homogeneous functional M x (-, z) in ( I16ab - 

The discrete analog of ( f6ab is by summation (for k = 1, 77 t) of the optimality con- 
ditions for (I16ab written at u — u k , i.e. 

f C4:e(v-4)dx+ f z\~ l Ku k T \v-u k T )dS > {f T ,v-u k T ) (18) 
Jn Jr c 

with e* from (I14ab and f T - f(£r), and tested by an arbitrary test-function v = v k . By 
comparison of values of (I16bl i at and an arbitrary 2, we get ${kr, u k , z k ) + ^(z^-z^ 1 ) < 
S{kT,u k ,z) + &Q(z-z k ~ 1 )- By the degree-1 homogeneity and the convexity of ^o(0, we 
further get the triangle inequality Moiz-z^ 1 ) < 3#o{z k T -z k ~ l ) + &o(z-Z k )- Re-organizing 
the first estimate and merging it with the second one, we obtain the discrete analog of the 
semistability i6b\ . namely: 

S{kT, u k T ,z k ) < S{kT, u k T ,z) + ^ (z-z k T -') - @o(4-4~ l ) ^ SQct, u k ,z) + @o(z-4), (19) 

A discrete analog of (l6cb as an inequality "<" can be obtained by testing the optimality 
conditions for d!6ab and d!6bb respectively by u k T -u\7 x and z k —z k ~ l (which, in fact, means 
plugging v = u\r l into (fT8l for the former test), and by adding it, benefiting from the 
cancellation of the terms ±$(kr, u k ,z k ~^) and by the separate convexity of §(t, •, •)■ This 
gives the estimate 

S(kr,uU) + t£^(^-, i^-) < W,«o,zo) + rj] (&£L,^). (20) 
i=i T T i=i T 

Let us by u XJ denote the continuous piecewise affine interpolant of the values {u k ) T k ^J Q , 
and by u x jT the piecewise constant "backward" interpolant, while T the piecewise con- 
stant "forward" interpolant. Analogously, we introduce z x , T and z interpolating values 





(z£)I5j> an d a l so fr and f T interpolating values {f T ) T Jj Q . In terms of these interpolants, we can 
write ( fT8l i. ( fT9l l. and (l20t more "compactly" as 



I Ce{xu x , T +u XtT ):e(v—u XtT ) dxdt + I z &^ T '(v-k^ t ) dSdf > I v— M^-^d? (21a) 
for any v € L 2 (I; H l (Q.; R d )) with v| 2c -n > 0, and 

U XiT (t), Z XlT (t)) < U x , T {t), Z) + ^o(H,tW). ( 2 lt>) 

for any z € L M (F C ) with < z < z x , T (t) on F c , and 

I ( I xCe ^,Ty- e (u x ,r) dx - { f T , u xT )^j dt + J a{za-z x , T (t)) dS 

+ £{t, u x , T (t),z x , T (t)) ~ AO, mo,Zo) =: ^,t(0 < 0, (21c) 

for any t = kr, k = 1, 77 t. 

Existence of (u\, z T ) solving (fl4l i is simply by a direct method applied to the underlined 
variational problems ([T6V Fixing ^ > 0, we can investigate the convergence for r — > 0. 
By a-priori estimates we have at disposal from (12 let , using Banach's and Helly's selection 
principles, we have immediately: 

Lemma 1. Assuming ( TTOt . (% T , Z^, T ) constructed recursively by (TT41 exists and, for^- > 
fixed, there is a subsequence (indexed by t's converging to 0) and u x e H l (I; H l (Q; M. d )) 
and z x 6 L°°(E C ) n V(7; Meas(r c )) such that 

u x ,r^u x mH\l;H l (Cl;m. d )), (22a) 

Z^t ^ z x in L°°(S C ) n BV(I; Meas(r c )), and (22b) 

z x , T (t) A z x (t) in L°°(r c ) for any t e [0, 7*]. (22c) 

Any (u x , z x ) obtained by this way is an energetic solution to (fl~|i due to DefinitionQ] 

In fact, the last claim required a limit passage in (|2H and then, a-posteriori, the proof 
of energy equality, which is rather technical and for details we refer to 11271 12911 . 

For further numerical study, cf . also Figures 0and |3]below, the important feature is that 
the residuum (£ XtT e L°°(7) in the discrete energy (im)balance (12 let can be controlled by 
making the time step r > sufficiently small: 

Proposition 2. Assuming again ( TTOt and x > fixed, it holds (even without any need of 
selection subsequences as in Lemma[T|i that 

lim lie^lbm = for any 1 < p < oo. (23) 

r— >0 

Moreover, for a selected subsequence satisfying (1221 . the weak convergence d22ab is, in 
fact, strong, i.e. in particular 

e(u XtT ) ^> e(ii x ) in L 2 (Q;R dxd ). (24) 

Proof. Essentially, for d23l . the only important point is to prove (124-b . Using ( 12 let for 
t — T, this can be seen from 

YCe(u x ):e(u x ) dxdt < liminf I ^Ce(M^ T ):e(M^ T )dxdf < lim sup I %Ce(u x T ):e(u x T ) dxdt 

Jq r->o Jq 

< AO,«o,zo) + limsup(^<f T ,^ T >df- £a(z -z x , T (T))dS - g(T,u x , T (T),z x , T (T))) 

<^(0,«o,zo)+ f <f,u*>dt- f a(zo-%(r))dS -*(r, ii^T), Z^Tj) 
Jo Jr c 

e(M^.):e(M^.) dxdf . (25) 



Note that, by d22b1 >. we have at disposal u x , T {T) -» u x {T) in //'(^; R d ) if X > is nxed > 
which we used together with (I22cl > for f = r to estimate lim sup T ^ —&(T, u XJ (T), z x , T (T)) < 
-£{J, u x (T), z x (T)) in d25b . Eventually, the last equality in (l25T > can be proved by limiting 
a regularization of the Signorini condition, cf. 11261 Step 4 in Sects. 8-9]. As a result, (l25l l 
proves 

lim I / ^Ce(il A , T ):e(i< A - T )djfdf = I ^Ce(M^):e(M^) dxdf. (26) 

T - >0 Je Je 

Using uniform convexity of the space L 2 (g; K^ x ' / ) equipped with the norm ||e|| := (J^ Ce:e dxdf) 1 ^ 2 , 
Za|> with d26l i allows for improvement of (122a) to the strong convergence 



u XtT ^>u x in H\l\ H l (Sl;M. d )), (27) 

hence also d24"l) is proved. 

Therefore, L j xCe(u XtT ):e(u XT ) dxdt in (12 lc) converges to f Q J n xCe(u x ):e(u x ) dxdt 
even uniformly in t. From d27i i, we have certainly u x l {f) — > M^(f) in H l (O; K c ') for any f, and 
using also ( 122c) . we have lim T ^o S(t , u XiT (f), z x , T (t)) = S{t,u x {f),z x {t)) for any t e [0, T]. 
Thus we have the convergence in ( I21cb for any f e [0, T]. As the sequence {€ x , T } T> o does 
not alternate sign and is bounded in L°°(7), by Lebesgue theorem it converges to some 
<£ x e L°°{J) strongly in LP (J) for any 1 < p < oo. Thus we showed that, in the limit for 

T — > 0, 

g(t,u x (t),z x (t)) + f f x^eiu^y.eiu^,) dx - (f,u x )dt- <f(0,wo,Zo) + f a(zo-Zy(0) d5 = ^(f)- 
Jo Jn Jr c 

(28) 

We used already in d25l l that the left-hand side of (|28| | is zero, i.e. here also € x = 0, so that 
(l23T l is proved. 

Eventually, we can realize that, in contrast to (l22t and (l24l i, the convergence d23l holds 
even for the whole sequence (indexed by a-priori chosen countable number of t's), which 
can be seen by a standard (contradiction) arguments based on uniqueness of the limit (here 
just 0). □ 

It should be remarked that, however, we did not prove validity of (|23T > for p — oo. 
Anyhow, even a coarser mode of convergence (f23) can ensure that the inequality ( 121 cb 
yields eventually the energy equality (113cl i, as we will pursue in what follows. 

Having the time-discrete viscous scheme, one may think about a convergence for both 
X — > and r — > simultaneously to obtain the Kelvin- Voigt-approximable solution to 
the quasistatic rate-independent problem. Here a certain circumspection has to be taken: 
obviously, lim^oA' c C<?(M r , T ):e(M A%T ) = for any t > fixed, cf. also Remark Q] below; in 
fact, this convergence is even strong in W l '°°(I; L l (D)). Therefore clearly, 

limlim^Ce(il^ T ):e(M A , T ) = 

and the energy balance ( 1 1 3ct would be obtained with p — as an inequality only; cf. also 
Figure [3] below. It is thus obvious that only some conditional convergence will lead to 
the desired p and the energy equality d 1 3cb as before. Obviously, we must consider rather 
lim^o lim T ^o- Linking (|27| | with d!2cb , we have 

w*-limlimx l Ce(u x T ):e(ii x T ) = w*-limx l Ce(u x ):e{u x ) = p, (29) 

^0 T ->0 x^O 



meant in Meas(g), and by (1231 we have obviously also 



lim lim T — 

X ->0t^0 x ' 



(30) 



meant in L P (I), 1 < p < oo. These two double limits can be merged under an implicit 
stability criterion & : R + -> R + so that both 

w*-lim xCe(u XtT ):e(u XtT ) = fx and lim <B XT — 0; (31) 

t<3T(x) T<.9( X ) 
r-»0,*-»0 T-M.x-tO 

cf. the arguments in the proof of [Ql Cor. 4.8(ii)]. In this way, we obtain a Kelvin- Voigt 
approximable solution according Definition [2] However, the stability criterion t < 3?(x) is 
not explicit and thus not of a direct usage in general. 



4. NUMERICAL IMPLEMENTATION AND COMPUTATIONAL EXPERIMENTS 

A general observation is that (fTol represents two recursive alternating linear-quadratic 
minimization problems which, after another spatial discretisation leads to linear-quadratic 
programming. On top of it, as no gradient of z is involved in S, ( 116bl i has a local char- 
acter and allows, after a suitable discretisation of F c , decoupling on particular boundary 
elements. Therefore, conceptually the proposed scheme leads to a very efficient numerical 
strategy for fixed x > an d t > 0. 

The essential difficulty is realization of the convergence (TJTl ). As the defect measure p 
is typically not known (and, on top of it, is not unique), we can hardly control the former 
convergence in OTb . Yet, we can at least control the latter one. To this goal, we devise the 
following conceptual algorithm relying on the convergence ( 1231 ) considered with p = 1 : 



(1) Set;f = xo > and t — tq > 0, and choose 
y > fixed. 

(2) Compute (u XT ,z x , T ) and ||^, r || L i (/ ). 

(3) If not ||(£^,tIIli(/) < Cx y , then put r := r/2 
and go to (2). 

(4) Put^ := x/2 and r := r/2. 

(5) If not^ < XtiBd, go to (2). 

(6) The end. 

Table 2. Conceptual strategy to converge with the 
viscosity x and the time step t to approximate the 
correct energy balance. 

In this way, we have at least the energetics in the limit under control if Xfimi would be 
pushed to zero. 

Proposition 3. The procedure from Table 2 is an algorithm in the sense that, for the pa- 
rameters xo > X&nai > 0, To > 0, C, and y > given, it ends after a finite number of loops, 
giving a solution with a viscosity^ smaller than the a-priori chosen xsmi- 

Proof. The only notable point is that, after finite number of refinement of the time discreti- 
sation, the condition ||Cy, T |lL 1 m ^ Cx y m Step (3) can be fulfilled. This follows from (l23l 
and the fact, proved in Proposition [2] that this holds for the whole sequence of the time 
partitions. □ 

4.1. Spatial discretisation by BEM 

To launch computational experiments, one naturally needs to perform still a spatial dis- 
cretisation. As <?(f, -, z) and ffl x (-, z) are quadratic functionals, (I16ab is a quadratic problem 
with the only constraint on F c which are linear. This allocation of all nonlinear effects 
exclusively on the boundary F c allows for using efficiently the boundary-element method 
(BEM) combined with linear-quadratic programming treating the variables on T c . 

BEM standardly uses so-called Poincare-Steklov operators which are known in specific 
static cases, here in particular for the homogeneous isotropic elastic material which we 



consider in what follows. Yet, we have to calculate the visco-elastic modification and here 
we benefit from choosing the ansatz of the tensor of viscous moduli as simply proportional 
to the elastic moduli, i.e. x<C. Therefore we can use BEM with the same Poincare-Steklov 
operators as in the static case only for a new variable v k := u\+x(u*—u*~ )/t; then, in terms 
of this new variable, one obviously has the Kelvin- Voigt strain € k = e(v k ), the velocity 



K- 



,,k-\ 



)/T 



(v*- 



,k-\ 



)/(t+x), and the displacement u k T = (tv£+^m* 1 )I(t+x)' which is to 



be used in (fT4l) . leading to the problem 

divCe<y r ) + f* =0 



< = 
t(e(v*)) 



t t (e(4)) + 4" 



T 



on £2, 
onr D , 

on r N , 



(32a) 
(32b) 
(32c) 



, TV k T +XU k T ' 
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-(' 



.TV k +XU k 1 



t n ( e (v A T ))+z*- 1 (l 



, TV T +XU K T 
T+X 



«)«) = 0, 

/ 1 

■)■/? > 0, 



(U g (v{)) + ^ 1 (K ^ + f^ 1 ).ir)((TV*+^- 1 )->0 = o, 



4 * 4" 



2(T +Af ) 2 



o*<ff, (4-4 _1 x^-ff) = o, 
(tv?+^- 1 )-(t4+^- 1 ) + %, u (4) 



on r r 



(32d) 



(tv* 1 +xu k 2 )/(t+x) proceeding recursively for k 



1, ...T/t € N. Then, like 



with m* -1 

(fTol l, one constructs the corresponding minimization problems in terms of {v k ,z k ). Also, 
evaluation of the energy balance (12 let in terms of v is possible at least approximately. More 
specifically, by using the Poincare-Steklov operator for the auxiliary variable v k which gives 
the equilibrium stress (in contrast to u k ), we calculate the test of the traction stress t(e(v k )) 
by velocity, i.e. the boundary integral 



-L 



,£-1 



,(t-i 



«*):«(«*) 



—Ce( U k T - l y.e(u k - 1 ) 



(33) 



where again T := <9Q, we obtain approximately the rate of stored energy and dissipation 
together, which can be used to express the overall energy balance as in ( 12 let at least ap- 
proximately by using boundary values only except the bulk contribution of / from (|5}, 
namely 

Ke(v x , T ))-ii x ,T dS - { f\, u_ XT )j df + J ^z x , T (t)Ku XiT (tyu XyT (t) - ^zqKu -u + a(zo~Z x , T (t)) dS 

- (frif), u XsT (t)) + <f T (0), u ) = £ x At) + t f f Ce(ux, T):e(u XtT ) dxdt, 

Jo Jo. 

(34) 

cf. d!7al > and (I21cb . The coefficient ^+r/2 in Q?b makes this expression only an esti- 
mate of the actual energetics (12 let which would need rather X- This additional term 
jT<Ce{ii X j):e{u XJ ) will vanish if r — > 0, being of the order 0(t/x)- Yet, it may not be 
entirely negligible for small x, which is a certain drawback of the BEM implementation. 

BEM also allows for avoiding transformation (|4j of the Dirichlet condition if u (or here 
rather v) is considered only on F c which is, due to (0, far from F D and thus the partial 
derivative <?/(■, u, z) have a good sense. In fact, the calculations presented below have been 
obtained by BEM implemented by a so-called collocation method. 



In what follows, we use this implementation in a two-dimensional geometry and, as 
already mentioned, isotropic material. In this situation, the Poincare-Steklov operator in- 
volved in BEM is well known; cf. Jl4l Sect.2.2]. As for the material, more specifically we 
use 

£ = 70GPa (Young modulus), (35a) 



JO (inSect.S2), , D . .. , 

V= \0.35 (inSect.E3, (Poisson ratto); (35b) 

thus, with 5jj standing for the Kronecker symbol, the elastic moduli tensor used in the 
previous sections takes the form 

vE E 
Qjkl = (l + v)(l-2v) tfytf " + 2^ i5ik6jl+ 6il6jk) - 

The viscosity of material describe by the relaxation time x will be varied and adjusted in 
particular cases below. 

4.2. Computational experiments: a simple test geometry 

In this section, we test the two-dimensional algorithm on a O-dimensional example 



from H29H where the viscous solutions as well as the limit for^ — » are explicitly known, 
together with a resulting nontrivial defect measure fi. We choose a rectangular specimen 
glued on one side and pulled on the opposite one by gradually increasing Dirichlet load 
in the normal direction, cf. Figure Q] Choosing the Poisson ratio makes the quasistatic 
problem essential O-dimensional (i.e. the strain, stress, dissipation rates, and delamination 
z are spatially constant). Such sort of tests are common in building geophysical models 
where it is called a one-degree-of-freedom slider. 

(. L = 1 00 mm h 

rigid obstacle 

\ I r N 



r, 



visco-elastic body 



\ adhesive T N 

Fig.fTJ Essentially a O-dimensional experiment (if the material is incompressible) with 
gradually increasing Dirichlet load and explicitly known solution. 

Considering still an isotropic adhesive K,j = Kdij, the initial condition uq = and zo = 1, 
and gradually increasing Dirichlet load w D (f) = v D t, we know analytically the solution 
of the viscous problem as well as the Kelvin- Voigt approximable solution of the inviscid 
problem (which is probably even unique). More specifically, there is a time, let us denote 
it by t RWtX or (for the limit \ — * 0) by f RUP , when the spontaneous and complete rupture 
happens; f RUPvr is determined only rather implicitly as a solution t of the transcendental 
equation (ao? + b x (l-s~ t ' t >')) 2 = 2a /K with the coefficients 

E , LEK E 

an — v n , by — —y rv n , and t y = y , (36) 

E+LK x (E+LK) 2 x E+LK ' 

and for ^ — * it converges monotonically to some limit, let us denote it by f RUP ; more 
specifically, 



'RUP,^ / RUP 




St mr = =-a/-=t. (37) 



For^ > 0, the response is given by 

1 for t < f, 



10 for t > ? RUPvr , 



and the stress cr x (t, •) = <r x (t) and the viscous dissipation rate xEe(u x ):e(ii x ) are constant in 
space (with values denoted by o~ x (f) and r x (t), respectively), while the displacement u x (t, •) 
is affine with u x (t, L) = w D (t) and u x (t, 0) = wM) with 

_ Uv D -a )t - bxil-e-ti'x) for t < t, 
w xW-\, ^-«-t^„)lx 



RUPv^ ' 



( w rup,AT v nA;up,#) e 



for ? > /„, 



^(0 = 



A = yc, h c 

L L 



(E+LK) 2 
X v 



vr 



-v D f, 



RUP 



L 

and w„ 



,- 2 ( f - f RUP^)^ 



RUP ' 
Vpf-Wy(f) 

Z ' 

for f < f, 
for t > t, 



RUPvV ' 



(38b) 
(38c) 

(38d) 



RUPvV ' 



where ao, b x , and t x are from (136) and w RUpa , := (v D -ao)f RU p^-^(l-e ' RUP,,r From d38db . 
one can see that the viscous dissipation rate^Ce(M Ar ):e(« A ,) concentrates in time when^ — > 
and, referring to (112ch the resulting defect measure ji takes the form 



u = (6, ® I) with S m 

H meas(Q) V w 



a — 

E 



(39) 



with 5, € Meas(Z) denoting the Dirac measure supported at t and 1 e Meas(Q) is the spatial 
constant measure with density 1 (i.e. the Lebesgue measure) on £1, and <? RUP is the energy 
stored in the bulk at the time of rupture t p „ when also the driving force d x = iKvtyvv^ 



reaches the activation threshold a; cf. 12911 for details about this calculation. Although yu 
is known and thus we could design the strategy from Table 2 to control also the difference 
xEe{ii XT ):e{u XJ ) - fi in some norm on Meas(0 which would be weakly* continuous, we 
intentionally do not want it because, in general (as also e.g. in Sect. l4.3l below). fi is not 
known. 

In addition to d35l l, we consider K - 150 GPa/m, a = 375 J/m 2 , v D = 267 yum/s, and T = 
0, 375 s. The length of the specimen is L = 0.1 m, as already depicted on Figure [1] while 
its cross-section is not important in this experiment. It is important that the implementation 
is able to hold the energetics with a good accuracy that can be efficiently controlled by 
making the time step small, as proved theoretically in Proposition|2]and shown on Figure|2] 
for a moderate selected viscosity x- The BEM spatial discretisation was coarse as all the 
quantities are either constant or affine in space in this "1 -dimensional" example, so the 
coarseness of the spatial discretisation is irrelevant. 
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Fig-El Illust ration of the time-dependent residuum -<E XiT (-) in the energy balance 
d2 let for r gradually decreasing as depicted from up to down, while x = 
6.25 x 10~ 3 s is fixed. The numerical error occurs especially around sudden 
rupture but is shown to converge to for r — > 0, as also proved in ((23}. 

The interplay between x and t and its influence on the energy balance is depicted on 
Fig-El clearly showing a very slow (resp. no) convergence for small^ > (resp. for^- = 0). 
The strategy from Table 2 chooses, in fact, a path decaying sufficiently slow from the left- 
upper corner towards the right-down corner in Fig.|3jleft): 
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Fig.|3] Left: 

the convergence of L'-norm of <£ X , T parametrized by x< documenting the theoretical result 
from Proposition|2]for p = 1. 

Right: 

L™-norm converges similarly in this example although this convergence is not theoreti- 
cally supported by Proposition^ 



For gradually vanishing viscosity Figure [4] displays respectively w x and <t x from (I38bb 
and d3=Scl i calculated numerically by a sufficiently small time step t. 
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Fig-El The strain (left) and stress (right) response; due to the symmetry, these tensors have only 
one nonzero component. 



Here, the defect measure /u is known from ( I39D ; now with f RLlp = 0.322 s and <? RUP = 
803.75 J/m 3 ; cf. d37l i and d39t . We can thus check the former convergence (|3T1 at least a- 
posteriori, which allowed us at least to tune the parameters for the algorithm from Table 2. 
Figure [5}left displays r x from d38dl i calculated numerically by a sufficiently small time 
step t. To visualize the weak* convergence to a Dirac measure, we display rather the 
overall energy dissipated by viscosity on the interval [0, t], i.e. L xCe(M A - iT ):e(«^ >T ) df, which 

should converge to J' /j. dt being just a jump at time f RUP of the magnitude S Rvr , cf . Figure[5}- 

right; again realize that spatial dependence is not interesting here as all these quantities are 

constant in space. § 

I 
F 



*=3 125n1s 

-#=t£? 

- X= 25 l 





3.0 

>. 2.5 

a 2.0 

I 1-5 
■I 10 
f 05 
1 °. 

■jj time 

Fig.|U Left: 

Convergence of the viscous dissipation rate r XtT = xCe(u XiT ):e(u XtT ) towards the defect mea- 
sure fi from £39}, i.e. here the Dirac at f RUP = 0.322 s for^ = 0.025x2"* with k = 0, 1,2, 3 
and decreasing r chosen according the strategy from Table 2, zoomed in and depicted on a 
selected time subinterval [0.3 , 0.375]. 

Right: 

Energy dissipated by viscosity over [0, f], i.e. JT xCe(u XtT ):e(u XiT ) dt, converging to 
the jump at f RUP = 0.322 s of the magnitude S RW = 803.75 J. Also the convergence 
t RVPX y / RUP from d37t is well documented. 



Remark 1. (Direct calculation of inviscid problem.) In principle, our semi-implicit time 
discretisation works for x — 0, too. Even, solving directly the inviscid problem is algo- 
rithmically much simpler. Yet, as pointed out at the end of Section [3] we cannot expect 
reasonable results if x will converge to too fast with respect to r, and in particular if 
straight ^ = would be used. Here, we saw it already on Figure [3] where, for^f = 0, the 
error in the energy balance practically remains constant no matter how the time discretisa- 
tion refines. On Figure [5] x — would cause all curves to degenerate simply to the f-axis, 
which shows fatal non-convergence of the overall viscous dissipation. Thus also the en- 
ergy balance cannot hold. It is surprising that u-, <x- and z-responses may still numerically 
converge to the correct solutions, as documented on Figure |6] which may seem to give a 
very efficient numerical strategy. We observed this phenomenon in all our calculations, in 
particular also on Figs.[T0land[T3lbelow; cf. also the discussion in Sect. [5] 
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Fig.[6] A comparison of the strain (left) and stress (right) response of a energetically 
justified small-viscosity solution with an unphysical result without any viscosity 
obtained by a semi-implicit formula; strongly zoomed in and depicted on a se- 
lected short time subinterval around rupture [0.320 , 0.324]: a surprisingly good 
match is achieved although energy does not match at all (since ji = without 
viscosity), cf. also Fig.f3]for^ = 0. 

Remark 2. (Stress- versus energy-driven rupture.) The rupture of Kelvin- Voigt approx- 
imable solutions is essentially stress driven, while the energy-driven rupture (i.e. energy 
dissipated by delamination is compensated by the elastic energy got from the bulk and 
adhesive, being related to so-called Levitas' maximum realizability principle and leading 
to so-called energetic solutions, cf. 11221 12311 ) occurs in general earlier, here in this simple 
example it would be at time ~\j2a(LK+E)/(\<lKE), as noted already in |29Tl . i.e. already 
at time 0.292 s. The stress-driven rupture seems to be much more natural (especially if a 
large bulk would lead to extremely early delamination) and is also preferred in engineering 
(where mostly the existence of solution and the calculations are not analytically justified, 
however), cf. Il7ll . or also the discussion about energy versus stress or global versus local 
minimization in mathematical literature 



4.3. Computational experiments: a fully 2-D example 

We now want to demonstrate applicability of the above developed methodology and 
algorithms to nontrivial situations where the defect measure /u is not known and typically 
is inhomogeneous, i.e. not distributed uniformly in space. Although we keep correct en- 
ergetics via tracking numerically the latter convergence in OTb . it should be emphasized 
that the calculations are not fully reliable because the former convergence in (|3~H cannot 
be checked. At this occasion, it should be however also emphasized that all the previous 
studies about defect measures have had been only purely theoretical and analytically mo- 
tivated (and being related, like here, with possible lack of regularity of weak solutions of 
various continuum-mechanical problems, exhibiting various concentration effects in con- 
trast to regular weak solutions where the defect measure vanishes, cf. J5, 12, 1210, 24]) 
and, except 12911 . existence of nontrivial defect measures had been rather only conjectured. 
In contrast to it, the simulations here represent, to our best knowledge, historically the very 
first attempt to see particular nontrivial spatially non-homogeneous defect measures. 

For our computational experiment, we use a similar geometry as in Section l4~2l but. in 
contrast to Figure Q] with a delaminating surface F c on a different side and (in our 2nd 
experiment) also different direction of loading, both intentionally breaking the symmetry 



considered previously in Section 14.21 We also consider more realistic Poisson ratio, cf. 
( l35l l. The speed of loading was taken the same in both experiments: w D (f) = v D t with |v D | = 
333.3 fim/s; the direction of v D was varied: horizontal in the 1st experiment and vertical in 
the second experiment, cf. Figure|7] In addition to (|33T >, we consider K = diag(X n , K t ) with 
K a = 150GPa/m and K t = 75GPa/m, and the fracture toughness a = 187.5 J/m 2 . We use 
X = 0.01s for all calculations 
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Geometry and boundary conditions of the 2-D problem 
considered. 



Calculations with t = 5xl0~ 3 s and 9.33xl0~ 3 s have been performed 300 time steps, up to 
T = 1,5 s and T — 2.8 s for the 1st and the 2nd experiment, respectively. Such T was big 
enough to achieve a complete delamination of the whole contact surface. 



4.3.1. Horizontal-loading experiment 

In contrast to the example from Section 14.21 we have now the traction force on T c 
nonhomogeneous, and it is interesting to see its evolution in time. This is depicted for 6 
selected snapshots, starting from t = 0.21 s with an equidistant step 0.025 s, on Fig.[8](upper 
and middle rows) which shows the delamination gradually propagating on F c from right to 
left. The displacement of the whole boundary has to be reconstructed by the Poincare- 
Steklov operators which is a conventional procedure in BEM and is depicted on Fig. [8] 
(lower row). 

f = 0.21 ( = 0.235 ( = 0.26 ( = 0.285 ( = 0.31 ( = 0.335 ( = 0.36 




Fig.ll Upper row: distribution of the tangential traction force in the adhesive along T c . 
Middle row: distribution of the normal traction force in the adhesive along Yq. 
Lower row: deformed configuration of gradually delaminating specimen under loading 
(1st experiment) from Fig.|7] the displacement depicted magnified 100 x horizontally and 
500 x vertically to make the vertical deformation more visible. 

To present spatial distribution of the defect measure, one must reconstruct the strain inside 
the domain Q.. This is a delicate (but anyhow executable) issue in BEM. To visualize the 
rate of viscous dissipation (which approximates the defect measure fi, cf. d!2cl >. and may 
exhibit time oscillations which would make visualization difficult), we display rather the 
overall dissipation on the interval [0, f], cf. Fig. [9] which approximates the total variation of 
the defect measure [//(-, x)](dt) as a function of x: 
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Fig.[9] The spatial distribution of the energy dissipated by (even very small) viscosity 
over the time interval [0, f], i.e. JC #Ce(w^ )T ):e(w^ [T ) dt depicted in a gray scale at 
6 selected time instances as also used on Fig.[8] 

It is not surprising that the dissipated viscous energy is bigger in the right-hand part of 
the specimen which is particularly stretched during the delamination. Perhaps noteworthy 
phenomenon is that this energy is not localized along the delaminating surface; we saw this 
effect already in the example in Section l431 where it was distributed over the whole volume 
uniformly. 

An analog of Figure Etright) displaying the force response t t~> JT t(e(w, u))(t, x) dS 
is on Figure [TOl left) together with a comparison with results obtained by the simplified 
inviscid algorithm from Remark [1] 
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Fig.fK)] Vertical and horizontal components of the reaction force on the Dirichlet loading for small 
viscosity^ = 0.01s and energy well preserved (left), compared with the inviscid solution 
calculated by semi-implicit method but energy balance completely violated as in Remarkf^ 
(right). As on Fig.[6] a surprisingly good match of this force response can be observed. 



4.3.2. Vertical-loading experiment 

Eventually, we briefly present most of the responses from Section 14.3.11 for another 
loading as indicated on Figure [7] Of course, the response is considerably different in 
some aspects, although the phenomena commented already for the horizontal loading are 
again observed. Now, the delamination propagates more slowly and we depict it with an 
equidistant step 0.45 s (instead of 0.025 s used in the 1st loading experiment) starting from 
t = 0.05 s: 



f = 0.05 t = 0.5 t = 0.95 t= 1.4 f = 1.85 t = 23 t = 2.75 




Upper row: distribution of the normal traction force in the adhesive along T c . 
Middle row: distribution of the tangential traction force in the adhesive along T c . 
Lower row: deformed configuration of gradually delaminating specimen under loading 
(2nd experiment) from Fig.[7] the displacement depicted 100 x magnified. 

The analog of Fig. [9] is on Fig.[T2j showing again that the viscous energy (and also the 
defect measure) can be supported in the bulk far away from the delaminating surface F c 
and here even a tendency to surprising symmetry in spite of nonsymmetry of the boundary 
conditions: 
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Fig-CH The spatial distribution of the energy dissipated by viscosity over 
[0,?], i.e. J xCe(u x , T ):e(u x , T ) dt depicted at 6 selected time in- 
stances as on Fig.f^ Surprising tendency to a symmetry even 
under nonsymmetry loading can be observed. 

Eventually, the force response corresponding to the previous Figures [TTI - [T2l is on Fig. fOT left) 
Like on Fig.|T0j there is again a surprisingly good match if calculated by the simplified al- 
gorithm from Remark[T]as depicted Fig.[l3jright), although there is no theoretical guaranty 
of such force-response match and obviously there is no match of energy. 
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Fig.HU Vertical and horizontal components of the reaction force on the Dirichlet loading (left) 
and its comparison with the simplified inviscid algorithm from Remark[T](right), again 
showing a surprising match as on Figures ©andf^ 



5. CONCLUSION 

We devised and tested a numerical strategy to approximate the natural notion of solu- 
tion to delamination of purely elastic material (as usually considered in engineering ap- 
plications). These solutions involve certain "defect measures" and follow the asymptotics 
arising from vanishing Kelvin- Voigt visco-elastic rheology which in the limit gives mere 
elastic material with the delamination driven naturally by stress rather than energy, as de- 



vised purely theoretically in 112911 without any time discretisation. 

We showed a delicate interaction between vanishing viscosity and time discretisation 
and difficulty to calculate physically relevant solutions. After testing the algorithm on a 0- 
dimensional example where exact solution is known, we calculated a couple of nontrivial 
2-dimensional examples by using BEM. Beside, we also compared the results with those 
obtained by a simplified inviscid algorithm ignoring defect measures and thus violating the 
energy balance, and showed a very good match of stress-strain responses in all investigated 
particular cases, cf. Fi gs. [6] [TOl and [13] Such an algorithm, called a Griffith model, was 
advocated already in 02511 although the desctruction of energy conservation was already 
pointed out there, and an investigation of some dynamical model leading to a correct limit 
in the quasistatic evolution advised, cf. lEM Sect. 3.2]. We conjecture that this simplified 
algorithm may converge to local solutions in the sense as introduced for a special crack 
problem in |31J and further generally investigated in IU9II . but the relation (and the phe- 
nomenon of good match) with the vanishing-viscosity approach remains still not justified. 

The importance of the above presented methodology for calculation such defect mea- 
sures would be pronounced in the full thermodynamical context like fll^, cf. also [0, 
Sect. 5.4], where the defect measure would naturally occur in the heat-transfer equation as 
a heat source and thus would influence temperature distribution inside the body and then 
backward the mechanics e.g. through thermal expansion or temperature dependence of me- 
chanical properties of the adhesive. Interesting observation from Figures [9] and [12] is that 
the defect measure (and, the possible heat production) may occur even in spots which are 
quite distant from the surface undergoing inelastic dissipative process of delamination and 
it is certainly difficult (or rather impossible) to guess its distribution by intuition. 



Acknowledgments: The authors thank Professor Alexander Mielke for discussion about the 
local-solution concept. 
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